Solving systems of linear equations on a quantum computer 
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Systems of linear equations are used to model a wide array of problems in all fields of science and engineer- 
ing. Recently, it has been shown that quantum computers could solve linear systems exponentially faster than 
classical computers 1 1 1, making for one of the most promising applications of quantum computation |2 |. Here, 
we demonstrate this quantum algorithm by implementing various instances on a photonic quantum computing 
architecture. Our implementation involves the application of two consecutive entangling gates on the same pair 
of polarisation-encoded qubits. We realize two separate controlled-NOT gates where the successful operation of 
the first gate is heralded by a measurement of two ancillary photons. Our work thus demonstrates the implemen- 
tation of a quantum algorithm with high practical significance as well as an important technological advance 
which brings us closer to a comprehensive control of photonic quantum information | 3 1. 



Systems of linear equations play an important role in var- 
ious fields, ranging from natural science and engineering to 
medicine and social science. The ability to solve such sys- 
tems underpins many modern technologies, including traffic 
flow analysis, computer tomography, and weather forecasting. 
Although systems of linear equations are an old problem — 
the two-dimensional version was investigated by the ancient 
Babylonians — it was only in 1811 that Gauss developed a 
general algorithm for solving them, Gaussian elimination 1 4l . 
Today, as the sizes of data sets are growing, Gaussian elim- 
ination is too slow and more advanced methods are needed. 
An increasing demand for a detailed understanding of larger 
and larger systems pushes the tools of classical computation to 
their limits. Modern data sets can be so enormous that finding 
a solution to a system of linear equations can be prohibitive 
even for the latest supercomputers. 

Quantum computers have attracted tremendous interest be- 
cause they could outperform classical computers at certain 
tasks [51. For a classical computer, the number of compu- 
tational steps needed to solve a linear system is at least pro- 
portional to the number of variables. By contrast, the new 
quantum algorithm could, in some cases, make the computa- 
tional time proportional to only the logarithm of the number 
of variables Q. An important difference is that the quantum 
algorithm calculates the expectation value of an operator as- 
sociated with the solution rather than the solution itself. Here, 
we demonstrate this algorithm using two full controlled-NOT 
(CNOT) gates acting on two qubits to determine the solution 
of a two-dimensional system of linear equations. We show 
various instances of the algorithm for systems of linear equa- 
tions with different characteristics. 



THEORY 

Solving a linear system of equations, given a matrix A and 
a vector 6, means finding the vector x such that Ax — b. If we 
rescale the vectors to ||6|| = ||x|| = 1, we can represent them 
as quantum states \b) and \x), and our initial task becomes 
finding the state \x) such that 



A\x) = \b). 



The solution we seek is 



\b) 

\\A-'m' 



(1) 



(2) 



We assume, without loss of generality |1|, that A is an 
N X N Hermitian matrix with eigenbasis and eigen- 

values {Aj}, and is rescaled so that < Aj < 1. The state \b) 
can be expanded in the eigenbasis, \b) = ^j^i and 
we aim to prepare, up to normalization, 



\x) 



^ 1 



(3) 



The algorithm consists of three main steps. Here, we de- 
scribe the basic idea of the algorithm; a more detailed discus- 
sion can be found in 1 1 1 and the Appendix). 

The first step is to apply phase estimation, a general proce- 
dure for decomposing quantum states in a particular basis 1 6 - 
[a. For this, we add an additional "eigenvalue register" of m 
qubits to our system, each initialized in the state |0). Phase 
estimation then transforms |6)|0)*^^ into Xl^i 
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FIG. 1: The simplest case of the quantum algorithm for solving systems of linear equations. Given a Hermitian matrix A and input outputs 
\x) — A~'^ \b) I \ A~^ I ^) II if the ancilla qubit is measured to be 1. (a) The complete circuit, as derived in the text, with JJ — exp 2^~^ A) 
and d — —2 arccos(Ai/A2), where Ai,2 are the eigenvalues of A and the integer n depends on the eigenvalues (see Appendix), (b) The local 

unitary R diagonalises A, A — ^ R. For the algorithm to work perfectly with this many qubits, Ai,2 must be such that JJ — 

R^ ^i) R = R^ ZR. This circuit reflects this simplification, (c) Optimised circuit. The middle qubit can be completely removed, and the 
controlled-rotation decomposed, to give the final circuit (see Appendix for details). \b) and R are arbitrary, (d) Experimental implementation 
of the circuit shown in c. The local unitary operations are implemented with the help of a combination of half- wave plates and quarter- wave 
plates. 



where the eigenvalues \Xj) are stored in the eigenvalue reg- 
ister to a precision of m binary digits. 

The second step is to implement the nonunitary map 
\Xj) Xj^\Xj). For this, we introduce an additional 



"ancilla qubit" initially in the state |1). Depending on the 
value of I Xj ) in the eigenvalue register, we implement a con- 
trolled Ry{Oj) rotation on the ancilla qubit. Here, Ry{0) = 
ex.-p{—i0ay/2) and ay denotes the usual Pauli matrix. With 
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FIG. 2: Experimental setup. Shown is the experimental implementation of two concatenated CNOT gates. The input is set by a polariser, 
which can be followed by a local unitary operation (LU). The two gates are connected by fibers. Different matrices A can be implemented by 
adapting the LUs, and different states \b) by adapting the input state. The figure shows the most general case of two concatenated CNOT gates, 
combined with general LUs. For the implementation of the algorithm, we chose some LUs to be the identity and obtained the case shown in 
Figure [TJi. 



Oj = —2 arccos(C/Aj), the state of our system becomes 




where C < miuj \Xj\. 

Finally, the third step is to run phase estimation in reverse 
to uncompute | A j ) , giving 



^/?,K.)(^^l-^|0) + -|l)j. (5) 

We measure the ancilla qubit, and if we observe a 1 we will 
have prepared \x) in the state register. If we know the eigen- 
values, we can maximise the success probability by choosing 
the largest possible C, C = miUj \Xj\. 

The runtime of the algorithm is O {\og{N) s'^ Hi'^ / e) , where 
s is the sparsity of the matrix, hz its condition number, and e 
the acceptable error 1 1 1. Furthermore, 0(/^^) can be reduced 
to 0{hi) with amplitude amplification ||9l. The best classical 
algorithms require 0{Ns^/l^log{l/e)) time, meaning that at 
constant s, tz, and e, the quantum algorithm is exponentially 
faster. On a quantum computer, determining all amplitudes 
in a quantum state scales exponentially with system size; in- 
stead, the strength of the algorithm lies in the determination 
of expectation values {x\M\x) of some operator M. 



The simplest case involves three qubits: one state qubit, one 
eigenvalue qubit, and the ancilla qubit. The state and eigen- 
value qubits replace larger registers which would be needed 
in a general implementation of the algorithm. Whereas using 
one state qubit means that 1 6) is a two- vector and A a 2 x 2 ma- 
trix, using one eigenvalue qubit imposes more subtle restric- 
tions. With one eigenvalue qubit, only a single binary digit of 
the eigenvalues is computed by the phase estimation, meaning 
that for the algorithm to work perfectly, it must be possible to 
distinguish the two eigenvalues with a single digit. Conse- 
quently, we choose the two eigenvalues to be of the form O.aO 
and O.al, where a is a sequence of binary digits. 

The complete circuit for the algorithm as described is given 
in Figure [T] which also outlines the procedure for optimis- 
ing the circuit to require only two qubits and two consecutive 
CNOTs acting on them (Figure [T]:). The circuit depends on 
the eigenvalues of A, the unitary R that diagonalizes it, 

•^ = ^*(oA°>. 

and the input state \b). The algorithm succeeds — i.e. the 
ancilla qubit is measured in the state |1) — with probability 
(Ai/A2)2. 

In our implementation, we choose different sets of eigen- 
values A = {Ai,A2}. The simplest case is a = 1, giving 
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the eigenvalues Ai = 0.10 = ^ and A2 = 0.11 = |. In 
this case phase estimation must read out the second bit of the 
eigenvalues. The procedure is analogous for the other sets of 
eigenvalues that we implemented (see Appendix). 



EXPERIMENTAL REALIZATION 

We implemented the algorithm using polarisation-encoded 
photonic qubits (Figures [T]l and [2]), where |0) and |1) denote 
horizontal and vertical polarisation, respectively fTO'.Tl |. Our 
implementation uses two different types of photonic CNOT 
gates lEMl. 

The first CNOT gate realizes the gate operation in a her- 
alded manner (191 EOl- It requires an entangled ancilla photon- 
pair and a measurement of two ancilla modes to herald that the 
gate has worked correctly. If two photons are registered in the 
ancilla modes, the gate has been successful without the need 
for a verification of the output state. Since the output photons 
need not to be measured, the application of a second CNOT 
gate is possible. 

Our second CNOT is implemented in a destructive way, 
where a coincident measurement of the output photons sig- 
nals the correct gate operation. The basic element of this de- 
structive CNOT gate is a polarisation-dependent beam splitter 
(PDBS) which has a different transmission coefficient T for 
horizontally polarised light (Th = 1) as for vertically po- 
larised light (Ty = 1/3) 1 14 1 . If two vertically-polarised pho- 
tons are reflected at this PDBS, they acquire a phase shift of 
TT. Two successive PDBSs with the opposite splitting ratios 
then equalize the output amplitudes. This setup, in combina- 
tion with two half- wave plates (HWPs) (see Figure [T]l) imple- 
ments a destructive CNOT gate. The gate has been successful 
if one photon is measured in each output mode. 

Combining these photonic CNOT gates with local unitary 
operations allows us to implement the circuit shown in Fig- 
ure [T] In our setup, we implement these local unitary oper- 
ations with the help of a combination of quarter-wave plates 
(QWPs) and HWPs (see Figures[T]l and[2]). A detailed descrip- 
tion of our experimental setup can be found in the Appendix. 



IMPLEMENTATION OF THE ALGORITHM 

We have implemented various instances of the algorithm, 
where we varied both the matrix A and the state \b). The 
matrix can be modified by tuning the local operation R (in 
which case the eigenvalues stay the same) or by adapting 0, 
which alters the eigenvalues of the matrix. The detection of 
the ancilla qubit in the state |1) announces a successful run 
of the algorithm and the preparation of the output qubit in the 
state \x). 

Figure |3] shows the results of a sample run of the algo- 
rithm. Two different matrices A were implemented by using 
different local operations Ri and R2. The eigenvalues were 
A = 1} in both cases. We then choose state vectors \bi) 
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FIG. 3: Experimental results, (a), (d) The figure shows two dif- 
ferent systems of linear equations, depicted by the matrices Ai 
and A2 as well as the state vectors \bi) and I62). (b), (e) The 
reconstructed density matrices of the experimentally obtained out- 
put state \x) are shown. These density matrices are obtained by 
choosing the local operations Ri = Rx{^7T).Ry(^7r) (b), and 
R2 = Rx{^7r).Ry{—^7T) (e). For both matrices, we choose the 
eigenvalues to be Ai = | and Ai = | by implementing Ry{0) as 
described in the main text. The fidelities of the reconstructed density 
matrices are 0.953 ±0.026 (b) and 0.976 ±0.010 (e). The wireframe 
shows the theoretical prediction, (c), (f) The quantum algorithm is 
based on determining the expectation value {x\M\x) of some opera- 
tor M with respect to the output state \x). Therefore, we also show 
the experimentally determined (blue) and theoretical (black) expec- 
tation values of several operators M. We choose the operator M to 
be the projection on the states |0), |±), and respectively, with 
1+) = (|0) + |1))/^, and |+.) = (|0> + i|l»/V2. 
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terns and to more precise and full phase estimation. In par- 
ticular, the ability to do arithmetic operations with quantum 
gates will enable the calculation of inverses (and thus the map 
\Xj) Xj^ \Xj)) on the fly. Even though current quan- 
tum computations in optical systems are proof-of-principle 
demonstrations, important insights can be obtained for future 
realisations in larger systems [ 2Tti23l . We are likewise hope- 
ful that our demonstration of this algorithm will enable future 
implementations of other, equally important algorithms that 
use it as a subroutine, including quantum algorithms for solv- 
ing nonlinear differential equations 1 24] and quantum data fit- 
ting 1251 . 



FIG. 4: The figure shows the solution of the system of linear equation 
for matrices with different eigenvalues. Experimentally, these are 
obtained by implementing different values of (see Appendix for 
details). For all matrices we run the algorithm for two input states 
\bi) = |1) = (0,1) and I62) = 1+) = (1,1)/V2. We achieve 
fidelities of 0.957 =b 0.010, 0.961 ± 0.013, 0.981 =b 0.009 for the 
input state \bi) (upper row from left to right) and 0.778 ± 0.031, 
0.773 ± 0.027, 0.832 ± 0.031 for the input state I62) (lower row, left 
to right). 

and 1 62) by setting different input states in our experiment (see 
Figure [3^). We run the algorithm and analyse the output state 
\x) via quantum state tomography. The density matrices of 
the resulting output states \x) are shown in Figure [s]:. As the 
algorithm relies on the determination of the expectation val- 
ues {x\M\x), the figure also shows expectation values with 
several operators M (see Figure [sji). 

Furthermore, we implemented the algorithm for a series 
of matrices A with three different sets of eigenvalues Ai = 
{^,f}, A2 = {|,|},andA3 = {f,|} (see Figure g. As 
explained in detail in the Appendix, the performance of the 
algorithm depends on the state R\b) that enters the gate. In 
order to analyse this behaviour, we chose two different input 
states, \bi) = |1) and I62) = |+) for each set of eigenval- 
ues, while keeping R equal to the identity matrix. The six 
resulting density matrices are shown in Figure [4] We achieve 
fidelities of up to 0.981 ±0.009 for and 0.832 ±0.031 for 
162)- This difference in fidelities arises due to the influence of 
higher-order emissions, which depends on the state R\b), as 
discussed in the Appendix. Additional data for a set of differ- 
ent input states and the three different choices of eigenvalues 
is shown in the Appendix. 



DISCUSSION 

The results presented here include the implementation of 
the simplest case of the quantum algorithm for solving sys- 
tems of linear equations, as well as the concatenation of two 
entangling gates acting on the same photonic qubits. We an- 
ticipate that increasing technological capabilities, including 
the implementation of more than two consecutive CNOTs, 
will allow the extension of the algorithm both to larger sys- 
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APPENDIX 
Theory 

Our work differs from the original proposal in 1 1 1 in sev- 
eral modifications that are needed to implement the algorithm 
with a limited number of qubits. Some of these modifications 
were mentioned in the main text, and we elaborate on them 
here, essentially describing the simplifications that are used to 
transform the circuit in Figure la to the circuit in Figure Ic. 

As described in the text, the most general circuit involves 
phase estimation, a controlled Ry rotation, and a reverse phase 
estimation (for a description of phase estimation in general, 
see Sec. 5.2 of [5]). If we restrict ourselves to one state qubit 
and one eigenvalue qubit, phase estimation is simply a con- 
trolled unitary between two Hadamard gates on the eigenvalue 
qubit (see Figure la). 

The unitary U that is chosen depends on what binary digit 
of the eigenvalue needs to be read out. Because we have as- 
sumed that A is scaled in such a way that its eigenvalues lie in 



the range (0, 1), we can express them as follows: 

Ai = O.xi, 1X1,2:^1,3^1,4 ... (7) 

A2 = O.X2, 1X2,2^2,3^2,4 . . . , (8) 

where the Xij are binary digits, or 1. If we wish to read out 

the n*^ digit, we must choose U = exp (27ri 2^~^A). Note 
its action on the eigenvectors: 

U\uj) = exp (27ri 2^-^^) \uj) (9) 

= exp (27ri 2^~"^0.Xj,iXj,2^j,3^j,4 . . .) (10) 

= exp (27rzO.X^-,n^j,n+l^j,n+2 . . .) \Uj) (11) 

^ exp (27rz 0.x j^n) (12) 

= {-lf--\u,). (13) 



The approximation of neglecting all digits past introduces 
errors to the procedure, which would be mitigated by using 
additional eigenvalue qubits. Without the additional qubits, 
we avoid the error by assuming, as in the main text, that the 
eigenvalues are of the form O.aO and O.al, differing only at 
the n*^ digit. 

The state and eigenvalue qubits, initialized to \b)s\0)E, 
transform thus: 

2 2 

|6>5|0)^; = E/5.>i)|0> ^ E/3.>.>T^ (14) 

2 

^^/3,>,)|%,). (16) 

The two qubits are now entangled, with each eigenstate in the 
state register accompanied by the n*^ digit of its eigenvalue. 

As described in the main text, = — 2 arccos(Ai/A2) in 
all cases. For example, if, as in Figure 3, Ai = 0.10 = ^ 
and A2 = 0.11 = |, in which case we read out the second 
digit, we get = —1.682. In the second column of Figure 4, 
Ai = 0.100 = ^ and A2 = 0.101 = |, and so i9 = -1.287 
(with the third digit read out). In the third column of Figure 4, 
Ai = 0.110 = ^ and A2 = 0.111 = |, and so 6> = -1.082. 

Once we assume that the eigenvalues differ in the n*^ digit, 
the action of U is seen to be very simple: it is merely a Z gate 
in the eigenbasis of A. That is, it leaves the first eigenstate 
unchanged and adds a phase of — 1 to the second. It can there- 
fore be written as /7 = R^ZR, where R diagonalizes A. It 
is with this fact, and remembering that a controlled- Z conju- 
gated with Hadamard gates gives a controlled-NOT, that one 
obtains the circuit in Figure lb. 

To get from Figure lb to Figure Ic, we first observe that the 
eigenvalue qubit can be removed. After the gate R, the first 
CNOT transfers the state of the state register to the eigenvalue 
register. This then operates the controlled-i?^(6>) rotation. In- 
stead of transferring the state first, it is equivalent to simply 
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(sin f, cos f) 
+> = (l,l)/v^ 
|C) = (cosf,sinf) 

|0) = (1,0) 



FIG. 5: The figure shows the different states R\b) which we have chosen as input states to our circuit. 



control the Ry{0) gate directly from the state qubit, giving 
this circuit: 



1^) -\R 
|i) — 



\x) 



Ry{e) 



Controlled single-qubit rotations have already been imple- 
mented in linear optics |[2T][22l[26l, but here we follow a dif- 
ferent approach. We decompose the controlled-i?^ rotation 
using the general method described in Sec. 4.3 of f5l, which 
immediately gives the circuit in Figure Ic. 



Experimental setup 

In our experiment, entangled photon pairs are produced by 
exploiting the emissions of a non-collinear type-II SPDC pro- 
cess 1271 . For this, a mode-locked Mira HP Ti:Sa oscillator 
is pumped by a Coherent Inc. Verdi V-10 laser. The pulsed- 
laser output (r = 200 fs, A = 789 nm, 76 MHz) is frequency- 
doubled using a 2 mm-thick lithium triborate (LBO) crystal, 
resulting in UV pulses of 0.75 W cw average. We achieve a 
stable source of UV pulses by translating the LBO to avoid 
optical damage to the anti-reflection coating of the crystal. 
The UV laser beam passes through a 2 mm-thick /3-barium 
borate crystal, gets reflected, and passes through the crystal a 
second time (see Figure [2]). The photons created during the 
first pass of the laser beam enter the first CNOT gate as the in- 
put (control and target) qubits. The state of these input qubits 
is modified using polarisers and additional local unitary gates 
that in principle allow for the creation of arbitrary input states. 
In our experiment, the control and target qubits are prepared 
in the states R\b) and respectively. Thus, we absorb the 
local operation R in the preparation of the input state. The 
photons created when the laser passes through the crystal the 
second time act as the entangled ancilla photon pairs which 
is required for the first CNOT gate. For this, we align our 
setup such that the entangled state |<l>+) = (|00) + \11))/V2 
is emitted ll20l . 

The photons interfere at the polarising beam splitters (PBS) 
as shown in Figure|4] The PBS on the control side (target side) 



is aligned such that it acts in the basis {|0), |1)} ({|+), |— )}). 
The photons are the filtered spatially and spectrally with the 
help of narrow-band filters (AA = 3nm) and by coupling 
them into single-mode fibers. A coincidence detection of the 
ancilla qubits in detectors 3 and 4 in the state | — )3 1 1)4 signals 
a successful gate operation. 

The output photons in the modes I and II (see Figure [4]) 
are then guided to the second CNOT gate. Wave plates be- 
fore and after the second CNOT implement the local rotations 
Ry{0/2),Ry{0/2),R,Sind R^ . 

Photons are coupled out, pass the polarisation-dependent 
beam splitters (PDBSs), and are coupled to single-mode fibers 
again. The success of the second CNOT operation is deter- 
mined by postselection on a coincidence detection in outputs 
1 and 2. The algorithm succeeds if the target photon is de- 
tected in state The output control qubit is in the state 
which is analysed by using HWPs, QWPs and polarising 
beam splitters; and a full state tomography of the output state 
\x) is performed. Errors are obtained from a Monte Carlo rou- 
tine assuming Poissonian counting statistics. These indicate a 
lower bound for the actual error that takes all the experimental 
imperfections into account. In our experiment, typical visibil- 
ities of the emitted Bell pairs are about 0.9 and higher-order 
emissions degrade the quality of our gate operations. The 
back-reflecting mirror is continuously moved back-and-forth 
to avoid any phase correlations between of the signal and the 
noise originating from higher-order photon emissions. Addi- 
tionally, imperfect visibilities on the order of 0.85 to 0.9 of 
the quantum interference at the PBSs in the first gate and the 
PDBS in the second gate contribute to errors. 



Experiment 

As mentioned in the main text, the fidelity of the output 
state \x) depends on the state R\b) that effectively enters the 
control input of the first CNOT gate. The Figures [6j [7] and [8] 
show the characterization of the state \x) for various input 
states. Figure [5] depicts the different input states we have cho- 
sen for our analysis. 

Our analysis shows that the fidelities of the obtained output 
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states vary from (64.7±4.2)% to (98.1 ±0.9)%. These varia- 
tions in fidelity arise due to the influence of higher-order emis- 
sions from spontaneous parametric down-conversion. These 
higher-order emissions can either occur when the beam passes 
the crystal the first time ("double-forward emission") or 
when the beam passes the crystal the second time ("double- 
backward emission"). 

If the input to the first CNOT gate is chosen to he R\b) = 
\0) or R\b) = \1) , 3. double-forward emission can never lead 
to a fourfold coincidence and thus signal a wrong "successful" 
operation of the first gate. In these cases, the photons entering 
the control input to the gate will either both be reflected or 
both be transmitted at the first PBS. For all other inputs to the 



gate, these higher-order emissions degrade the fidelity of the 
output state as demonstrated in our analysis. 

The "double-backward emission" can in principle never 
lead to a fourfold coincidence because the photons can never 
split up due to quantum interference. However, due to the vis- 
ibility of the entangled Bell pairs of 0.9, this quantum interfer- 
ence does not work perfectly. About 10% of our total counts 
arise from these events, which also influences the fidelity of 
the output state \x). However, for the input R\b) = |1) the 
fidelity does not seem to be affected. In this case, the noise 
cannot be distinguished from the signal. The influence of the 
double-backward emission increases for the other input states 
and reaches a maximum f or R\b) = |0). 




FIG. 6: The figure shows the different states R\b) which we have chosen as input states to our circuit. 




FIG. 7: The figure shows the different states R\b) which we have chosen as input states to our circuit. 




FIG. 8: The figure shows the different states R\b) which we have chosen as input states to our circuit. 



